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Abstract: We are interested in the modehng of wave propagation in an infinite bilayered 
acoustic/poroelastic media. We consider the biphasic Biot's model in the poroelastic layer. 
The first part is devoted to the calculation of analytical solution in two dimensions, thanks 
to Cagniard de Hoop method. In this second part we consider the 3D case. 
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Solution analytique pour la propagation d'ondes en milieu 
stratifie heterogene acoustique/poroelastique. Partie II : en 

dimension 3 

Resume : Nous nous interessons a la modelisation de la propagation d'ondes dans les mi- 
lieux infinis bicouche acoustique/poroelastique. Nous considerons le modele bi-phasique de 
Biot dans la couche poroelastique. La premiere partie est consacree au calcul de la solution 
analytique en dimension deux a I'aide de la technique de Cagniard-De Hoop. Dans cette 
deuxieme partie nous considerons le cas de la dimension 3. 

Mots-cles : Modele de Biot, ondes poroelastiques, ondes acoustiques, couplage acous- 
tique/poroelastique, solution analytique, technique de Cagniard de Hoop. 



Analytical solution 



Introduction 

The computation of analytical solutions for wave propagation problems is of high importance 
for the validation of numerical computational codes or for a better understanding of the 
reflexion/transmission properties of the media. Cagniard-de Hoop method [U |6] is a useful 
tool to obtain such solutions and permits to compute each type of waves (P wave, S wave, head 
wave...) independently. Although it was originally dedicated to the solution of elastodynamic 
wave propagation, it can be applied to any transient wave propagation problem in stratified 
media. However, as far as we know, few works have been dedicated to the application of this 
method to poroelastic medium. In [12] the analytical solution of poroelastic wave propagation 
in an homogeneous 2D medium is provided and in [13] the authors compute the analytical 
expression of the reflected wave at the interface between an acoustic and a poroelastic layer 
in two dimension but they do not explicit the expression of the transmitted waves. 

In order to validate computational codes of wave propagation in poroelastic media, we 
have implemented the codes Gar6more 2D [10] and Gar6more 3D [11] which provide the com- 
plete solution (reflected and transmitted waves) of the propagation of wave in stratified 2D 
or 3D media composed of acoustic/acoustic, acoustic/elastic, acoustic/poroelastic or poroe- 
lastic/poroelastic The 2D code and the 3D code are freely downloadable at 

[iittp : //www . spice-rtn . org/lib rary/sof tware/Gar6more2D[ 
and 

|http : //www. spice-rtn. org/library/sof tware/Gar6more3D[ 
In previous studies [HI [9] we have presented the 2D acoustic/poroelastic and poroelastic/poroelastic 
cases and we focus here on the 3D acoustic/poroelastic case,the 3D poroelastic/poroelastic 
case will be the object of forthcoming papers. We first present the model problem we want 
to solve and derive the Green problem from it (section 1). Then we present the analytical 
solution of wave propagation in a stratified 3D medium composed of an acoustic and a poroe- 
lastic layer (section 2) and we detail the computation of the solution (section 3). Finally we 
illustrate our results through numerical applications (section 4). 

1 The model problem 

We consider an infinite three dimensional medium (O = R ) composed of an homogeneous 
acoustic layer Q~^ = R x] — cx),0] and an homogeneous poroelastic layer Q~ = R x [0, -|-oo[ 
separated by an horizontal interface T (see Fig. [1]). We first describe the equations in the 
two layers ( §l.ll and §1.21 ) and the transmission conditions on the interface T ( §1.3|) . then we 
present the Green problem from which we compute the analytical solution ( §1.4p . 
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Figure 1: Configuration of the study 



1.1 The equation of acoustics 

In the acoustic layer we consider the second order formulation of the wave equation with a 
point source in space, a regular source function / in time and zero initial conditions: 

p+ _ v+^AP+ = 5^5y5,_hf{t), in f^+x]0,r], 



1 



in Jl+x]0,T], 
inSl+ 
_ U+(x, y, 0) = 0, U+(x, y, 0) = 0, in 0+ 



P+{x,y,0)=0,P+{x,y,0)=0, 



(1) 



where 



P"*" is the pressure; 



[/"*" is the displacement field; 



V~^ is the celerity of the wave; 



p'^ is the density of the fluid. 
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1.2 Biot's model 

In the second layer we consider the second order formulation of the poroelastic equations [U 



HE] 



p-tj^ +pjW - V ■ s- = 0, 



1 



inS7-x]0,T] 
in n-x]0,T] 



pju,+p;,w +^w + vp- = 0, 

S- = \-V -U- I^ + 2p-e{U-) - f3- P- I-i, inS7-x]0,T; 
1 



m 



P- + (3- V -U- +V -W- = Q, 



[/7(x,0) = 0, W-{x,Q) =0, 
L/7(x,0) = 0, VF"(x,0) =0, 



in S7-x]0,T], 



in il 
in O" 



(2) 



with 



3 QY,-. 
(V • S^)i = ^^ "o"^ Vi = 1, 3, /a is the usual identity matrix of A^3(]R), 

. C/X J 



and £{U g ) is the solid strain tensor defined by: 



e^J{U) 



7^ + 



2 \dxj dxi 



In ([2]), the unknowns are: 



• Ug the displacement field of solid particles; 

• W = (j)~{U^ — Uj), the relative displacement, C/7 being the displacement field of 
fiuid particle and 4>~ the porosity; 

• P^ , the fiuid pressure; 

• S^, the solid stress tensor. 

The parameters describing the physical properties of the medium are given by: 

• p^ = (p^ p^ + {1 — 4>^)p^ is the overall density of the saturated medium, with p~ the 
density of the solid and pT the density of the fiuid; 



Pyj = a Pf/cf) 1 where a the tortuosity of the solid matrix; 
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/C = n /rj , K is the permeability of the sohd matrix and r/ is the viscosity of the fluid; 



• m~ and /3~ are positive physical coefficients: /3~ = 1 — K^ jK^ 

r - 1 ^-'- 

and wT = cp^ /Kr + (/3~ — (t)^)/Kj , where K^ is the bulk modulus of the solid, 
KJ is the bulk modulus of the fluid and Kj^ is the frame bulk modulus; 

• n' is the frame shear modulus, and A^ = K^ — 2/i^/3 is the Lame constant. 

1.3 Transmission conditions 

Let n be the unitary normal vector of T outwardly directed to. The transmission conditions 
on the interface between the acoustic and porous medium are [5] : 

( W n = (U+ -U-)-n, 

P-=P+, (3) 

Y.~n = -P+n. 

1.4 The Green problem 

We won't compute directly the solution to ([1][2][3|) but the solution to the following Green 
problem: 

p+ -V+^/\p+ = 5^5y5,.h5u inl7+x]0,r], (4a) 

w+ = — ^Vp+, inO+x]0,r], (4b) 

p" it^ + pj ii>~ - V • cr~ = 0, ini7"x]0,T], (5a) 

pjw7 +p^^" + -:^w" + Vp" = 0, inf]"x]0,T], (5b) 

a- = \-V -u- I-i + 2^.~£{u-)- 13-p^ I:i, \n ^-y.]{),T], (5c) 

1 



m 



■p^+ZJ^V-it^ + V-w" =0, inJl"x]0,T], (5d) 

w ■ n = (w^ — uj) ■ n, on F, (6a) 

p^ = p^, on r, (6b) 

a~ n = —p^ n, on T. (6c) 

The solution to ([T]12]13I) is then computed from the solution of the Green Problem thanks to 
a convolution by the source function. For instance we have : 

t 
+ f 



P+{x,y,t)=p+{x,y,.)*f{.)= / p+{x,y,T)f{t-T)dT 

Jo 
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(we have similar relations for the other unknowns). We also suppose that the poroelastic 
medium is non dissipative, i.e the viscosity r]~ = 0. Using the equations ([5cl I5dp we can 
eliminate a~ and p" in Q and we obtain the the equivalent system: 



p" u~ + Pj ih' - a- V(V • uj) + p" V x (V x uj) - m~P~V{V ■ w) = 0, z < 



prU^+p^ib — m j3 V{y-Ug) — m ViV -w ) = 0, 



z< 



(7) 



_ _ _ _ _2 

with a" = A^ + 2p^ + m^ f3~ . 



And using the equation (I4b|) the transmission conditions ([6]) on z = are rewritten as: 



ilsz + w^ 



1 



dzP^, 



—m 13 V -Ug — m Vw = p~^ , 
dzU~^ + dxuj^ = 0, 
dzu-y + dyu-^ = 0, 

_ _2 _ _ _ _ 

(A~ + 7n~f3~ )V ■ m7 + '^IJ'dzUl.^ + m^ (3" V • w 



-P 



8b 
8c 
8d 



(9 



We split the displacement fields u^ and Ur into irrotationnal and isovolumic fields (P-wave 
and S-wave): 

w7 = ve- + vx*- ; i,;- = ve^ + vx*^. 

The vectors ^u and *^ are not uniquely defined since: 

vx (*f + VC) = vx *f, yee{u,w} 

for all scalar field C. To define a unique *£ we impose the gauge condition: 

V • *f = 
The vectorial space of ^^ verifying this last condition is written as: 



^f 



dy 

-dx 




^Ll + 



52 

^xz 

-52 -92 



^ 



^,2) 



where ^ i^i and ^^^2 are two scalar fields. The displacement fields u^ and w are written in 
the form: 



u. 



ve- + 



52 
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w 



ve- + 



5'. 




92 - 

'-'xx 


52 

yy 


dl. 




dl 




-dlx- 


dly 



^ 



uA 



5, 



y 

dr. 



^ 



wA 



A^ 



«,2 



(10) 
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We can then rewrite system ([7]) in the following form: 



A-Q- -B-Ae- 


= 0, 


z<0 


Ki - ^s'^K,i 


= 0, 


z <0 


^^,2 - Vs'a^-, 


= 0, 


z <0 


^ 1/1 ^ i; 1 




z <0 



(11) 



Pw 



where = (0^ , 0^)*, A and B are 2x2 symmetric matrices: 



A- 



P Pf 

Pj Pw 



and 



B- 



v; 






m 



PPw 



\ P Pw- Pf 



is the S-wave velocity. 



We multiply the first equation of the system (llip by the inverse of A. The matrix A^^B 
is diagonalizable: A^^B = VDV^^, where V is the change-of-coordinate matrix, D = 
diag{Vfr , Vf^ ) is the diagonal matrix similar to A~^B, Vfr and Vf^ are respectively the 
fast P-wave velocity and the slow P-wave velocity (Vps < Vpf)- 



Using the change of variables 

we obtain the uncoupled system on fast P-waves, slow P-waves and S-waves: 



(12) 



/■ -T- 



DA^- 



0, 



u, -V7 A*-, = 0, i = l,2, z<0 






Using the transmission conditions (|8cj) - ()8dp . we obtain: 



z <0 



(13) 



Pw 



z <0 

= 0, 
= 0, 



on r, 
on r, 



(14a) 
(14b) 



with A_|_ = d'^^ + dy . Applying the derivative dy to the equation p4a|) , dx to the equation 
(I14bp and subtracting the first obtained equation from the second one, we get: 



(a,Ax)A^-2 = 0' on r> 



(15) 
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moreover, using that ^„ 2 satisfies the wave equation: 



* 



u,2 



^S ^K,2 = 0, ^ < 



and that Ug and w satisfy, at t = 0, u^ = u^ = w = w = 0, we obtain: 



K,2 = 0, ^ < 0, 



and from (I14p we deduce the transmission condition equivalent to (|8cl) and flSdh : 

2d!e- + {dl-A^)^-, = o, onr. 

Finally, we obtain the Green problem equivalent to (|4l5|6p : 



(16) 



f p+-V+^Ap+ = 6,6y6,^hSt, 



z>0 



<t>r - Vr^A<i>r = 0, i£ {Pf, Ps, S} z<0 



(17) 



where we have set ^g = ^^ ^ in order to have similar notations for the Pf, Ps and S waves. 
The operator B represents the transmission conditions on T: 



B 






^a, (ru+r2i)d%, 



{Vi2 + V22)dlu (^-i)4Ax 

Pw 



. m-{(i'Vii+V2i) ^2 m-{(3-Vi2 + V22) r,2 

r- 2 tt ^r-2 tt 




1 



2Vu d. 



Vps 

2Vi2 d. 



B. 



'42 



B. 



43 



dt-A 



-2/i-a,Aj 



P 



Pf 



Ps 



S A 



where Vij, i,j = 1,2 are the components of the change-of-coordinates matrix V, B42 and B43 
are given by: 



7-2n 



B„ = (A-+.n-rW.+>..-rP.. g. ^ 2^_^^^3.^ 



V", 



_ 2 
Pf 



B. 



'43 



(A- + m-l3'^)Vi2 + m-p-V22 



V, 



dl + 2^x~Vl2^l,. 



Ps 



To obtain this operator we have used the transmission conditions ()8a|8b|16|8ep . the change 
of variables Q and the uncoupled system (fTBj) . 

Moreover, we can determine the solid displacement u^ by using the change of variables 
dni) and the fluid displacement it"*" by using (I4bp . 
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2 Expression of the analytical solution 

Since the problem is invariant by a rotation around the z-axis, we will only consider the case 
y = and x > 0, so that the y-component of all the displacements are zero. The solution for 
y 7^ or X < is deduced from the solution for y = by the relations 



p{x, y, z, t) = p{yJliP^Ty^, 0, z, t) 



Usxix,y,z,t) 



Usy{x,y,z,t) 



-s/oF+y^ 



\l x^ + y^ 



=u, 



,x(v^r2 + y2,0, 2;,i) 



,(Vx2 + y2,o,z,t) 



Uszix,y,z,t) = Usz{Vx^-HJ^,0,z,t) 
To state our results, we need the following notations and definitions: 



(18) 
(19) 

(20) 
(21) 



1. Definition of tiie complex square root. For q^ G C\IR , we use the following 
definition of the square root g{qx) = qj '■ 

9{Qxf = qx and ^e[g{qx)] > 0. 

The branch cut of g{qx) in the complex plane will thus be the half-line defined by 
{qx G ]R~} (see Fig. [5]). In the following, we'll use the abuse of notation g{qx) = i y/—Qx 
for qx £ IR . 



9m(g,)> 


^ 


+7T^ 


• 1 ^ 


-V 


J ^e{q 



Figure 2: Definition of the function x i-^ (x)^' ^ 



2. Definition of the fictitious velocities For a given (/ G IR, we define the fictitious 
velocities V~^{q) and V^{q) for i G {Pf,Ps,S} by 



V+ := V+(g) = V^ 



1 



1 + y+2fl2 



andVr:=V-{q)=Vr 



1 + Vr^q^ 
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These fictitious velocities will be helpful to turn the 3D-problem into the sum of 2D- 
problems indexed by the variable q. Note that V^(0) and Vj"(0) correspond to the real 
velocities V~^ and V^^ . 

3. Definition of the functions k+ and k^ . For i G {Pf,Ps,S} and {qx,qy) G C x K, 
we define the functions 



K+ :=K+{qx,qy) = f 77; 



and 



i^i ■= Hi {qx,qy) 



K- 



I 2 ' 1x ' 1y 



+ ql + "^ 



1/2 



1/2 



1 



2/ ^ + 'ix 



2 ' ^x I Hy 



V+\qy) 






1/2 



1/2 



4. Definition of the reflection and transmission coefficients. For a given {qx,qy) £ 
C X IR, we denote by U{qx,qy), Tpf{qx,qy), Tpsiq-x^Qy) and Ts{qx,qy) the solution of 
the linear system 



^{qx,qy) 



T^{qx,qy) 
'^pf{qx,qy) 

Tps{qx,qy) 
Tsiqx^qy) 

where the matrix A{qx-, qy) is defined by: 
'^{qx,qy) = 

i<^{qx,qy) 



2K+{qx,qy)V+^ 



i^'^{qx,qy) 

1 



1 



(22) 



Pf 



-.2 I „2\ 



1 



1 



(Pii + 'P2i)Kpj{qx, qy) {Vu + 'P22)Kps{qx, qy) 1 ] {qx + q-. 



m 



m 



V, 



{P~Vu + P21) —^ {P~Vu + 7^22) 



2VuKpj:{qx,qy) 

M,2iqx,qy) 



V, 



Ps 

2V12K p^{qx,qy) 

M,?,{qx,qy) 



i^s {qx,qy) + ql + ql 
2p~iql + ql)Ksiqx,qy) 



with 



A4,2{qx,qy) 



Ai,3{qx,qy) 



X- + m-p-^ Vu + m-(3-V2i 



y, 



_ 2 
Pf 



2 

+ 2fi^K,~j, {qx-,qy)V_ 



11) 



A" + m-(5-" P12 + m-/3-p22 



F, 



_2 
Ps 



_ 2 

+ 2^'Kp^ {qx,qy)Vl2. 



We also denote by T^nax the greatest velocity in the two media: Vmax = max(F+, Vp r, Vp^, V^. 
We can now present the expression of the solution to the Green Problem : 
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Theorem 2.1. The pressure and the displacement in the top medium are given by 



p+(x, 0, z, t) = P^^^{x, z, t)-\ — -^{x, z, t) and ti+(x, 0, z, t) = u^^^{x, z, r) dr+u^^^x, y, r) dr, 
and the displacement in the bottom medium is given by 

U~{x,0,Z,t) = Upj^{x,Z,t) + Up^{x,z,t) + Ug{x,Z,t) 

where 



• p- and u- are respectively the pressure and the displacement of the incident wave 
and satisfy : 



pi" (x,z,t) 



^inc.xi^'^'*^ •" 



^inc.z^'^' ^' ''^ •" 



S{t - to) 

ATrV+'^r 

xtH{t - to) 
47ry+V3p+ 

{z - h)tH{t - to) 



47ry+V3p+ 

where 5 and H respectively denote the usual Dirac and Heaviside distributions. More- 
over we setr = (x^ + (2; — /i)^)^'^ and to = r/V^ denotes the time arrival of the incident 
wave at point (x,0,z). 



* ^^ f and u'^ r are respectively the primitive of the pressure with respect to the time and 
the displacement of the reflected wave and satisfy : 



Ce/^' ^' *^ 



u'^ ,■ (x, z,t) 
refx^ ' ' ' 



U^ r ix, Z,t) 

refz^ ' ' ^ 



qi(t) Qm 



K+{v{t,q))n{v{t,q)) 



vrVvV+^IW 



dq, 



qi{t) 9m 



iv{t,q)K+{v{t,q))TZ{v{t,q)) 



qi{t) Qm 



K+\vit,q))nivit,q)) 



dq, 



Ti'^rp+yjq^ +q^{t) 



dq, 
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X V^ 
if th < t < to and — > 

"-1 r Vr„, 



^rej^X' ^' *) 



qi(t) Qm 



K+{vit,q))nivit,q)) 



IT- 



'r^/q'^ -qi{t) 



<?0 



(t) ^e 



K+{jit,q))ni^it,q)) 



vrVy^p) 



d(7 



dq, 



u^ ,■ (x, z, t) 



<ji(t) Qm 



(i^;(t,(?)/t+(^;(t,(?))7^(^;(t,g)) 



go(t) 



TT' 



V/9+Vg2_^2(^) 



dq 



+ 



90 



(t) JRe 



ij{t,q)K+{-f{t,q))TZ{-f{t,q)) 



U^ r (X, Z,t) 



qi{t) Qm 



TT'^rp+y/q'^{t) -g2 



K+'(i;(t,g))7e(i;(t,g)) 



dg, 



go(t) T^'^rp+yjq^-q^{t) 



dq 



+ 



go 



(t) Jfte 



/.+^(7(t))7^(7(t)) 



vrVp+vWiW 



dq, 



X V^ 
ifto<t< tu and - > —- 



Ce/^' ^' *^ 



It"*" <■ (x, z, t) 






go 



(t) Sfte 



«+(7(t,(?))7^(7(t,(?)) 



vr- 



vvi(t)^? 



dq, 



Qo 



(t) 'Sie 



■n{t,q)K+{jit,q))n{j{t,q)) 



2^^rp+y/q'^{t)-q^ 



dq. 



q»{t) Re 



K+\^{t))n{^{t)) 



vr- 



■rp+^/qQ{t) -q^ 



dq, 



X 



V+ 



if t^ < t and — > — — or ifto<t and — < 



V+ 



and 



f ''max ^ ''max 

^^gdx, y,t) = and u^^Ax, y,t) = else . 

We set here r = (x^ + (z+/i)^)"'^'^ andt^ = r/V^ denotes the arrival time of the reflected 
volume wave at point (x, 0, z), 



thi = (z + h)^ 



1 



V- 



1 |x| 

* rn a -V- ' TTl 



max * max 

denotes the arrival time of the reflected head-wave at point (x, 0, z) and 



thn 



1 



1 



z + hV y+2 T/2 



(23) 



(24) 
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denotes the time after which there is no longer head wave at point {x, 0, z), (contrary to 
the 2D case, this time does not coincide with the arrival time of the volume wave). We 
also define the functions 7, v, qo and qi by 



xt z ~\~ h / 1 
'y:{teIR\t>to}xm^C:='y{t,qy)=i— + 



r V r^ V+\qy) 



V : {t e IR\thi <t < t/i2}x^^ C := v{t,qy) = -i 

qo : m ^ m := qo{t) = 



z + h 



t2 



+ ^t 






V+\qy) r^ r 



and 



qi: m^ IR:=qiit) 



\ 




t- (z + h). 



1 1 



I/+2 V^ 

^ ' ma 



y2 • 



Remark 2.1. For the practical computation of the pressure, we won't have to explicitly 
compute the derivative of the function ^^ r (which would be rather tedious), since 

r,+ X: f = dtE^ * f = £^ * f' 
ref ^ref ^ref 

Therefore, we'll only have to compute the derivative of the source function f . 
• UpAx,z,t) is the displacement of the transmitted Pf wave and satisfies: 



Pii /"^i^*) 

Upf^{x,z,t) = 5- / 3f?e 

T^ Jo 



ifthi <t<to and \Qm[y{to,0)]\ < 



^e 
1 



•'ma: 



dv 
yv{t,q)Tpf{v{t,q)) — {t,q) 



dv 



dq, 



Kpfivit, q))Tpf{v{t, Q))-Q^it, q) 



dq, 



^ Jo 



^e 



Vn r^^*^ 



^e 



Upj^^{x,Z,t) 






-7r2 , 
^ ./O 



^e 



+ ^r"\e 



vr^ 



go(t) 



n(.t,q)Tpf{'y{t,q)) — (t,q) dq 

dv 
iv{t,q)Tpf{v{t,q)) — {t,q) dq, 

_ d^ 

i^Pfhit, QWpfhit, ^))^(*' «) 

_ dv 

Kpj{v{t, q))Tpf{v{t, Q))-Q^{i, q) 



dq 
dq, 
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iftQ<t<th^ and \Qm[y{to,0)]\ < 



yxr 



Upr (x,Z,t) 



^11 /"*'*•*-' 



-7r2 , 



3fie 



Upj,^{x,z,t) = ^TT I ^e 



-7r2 , 



[7(t,g)rp/(7(t,g))|^(t,(?) 



dq, 



i^'Pfiiit, qWpfiiit, 9))-^(*' 9) 



dg, 



i/t/i2 < t and |9m [7(^0, 0)]| < — or if to < t and \Qm[y{to,0)]\ > — 

''max ''max 

and Upr{x, z,t) =0 else. 

to denotes here the arrival time of the Pf volume wave at point (x, 0, z) (we recall in 
appendix the computation of to), 



thi = h\ 



1 



1 



1 



T/- ^ V^ ' K,«x 

Vpp "^max ►max 



y ' max 

denotes the arrival time of the Pf head wave at point (x,0,z), 



1 \x\ 

+ 



(25) 



thn 



h _z_ 

Cl C2 



(26) 



denotes the time after which there is no longer head wave at point (a;, 0, z), where 



1 



1 



Cl 



T/+2 y2 



and C2 



1 1 



T/2 T/2 • 

''p/ "^max 



The function qo : [to J +00] h^ iR^ is i/ie reciprocal function of to : iR^ ^^: [toj+oo], 
where to{q) is the arrival time at point (x, 0, z) of the fictitious Pf volume wave, prop- 
agating at a velocity V'^{q) in the top layer and at velocity VpAq) in the bottom layer 

(we recall in appendix the computation ofto{q)). 
The function qi : [ti ; to] ^^ ^^ is defined by 



idt) 



\ 



1 1 

t + Z 



V p r * ma 



T/+2 F2 
i' ^ ma 



y2 • 



The function 7 : {(t,^) G iR''' x IR^ \ t > io{q)} ^^ C is implicitly defined as the only 
root of the function 



^h,Q,t) = -z 



+ 7M +h( —IT— + 7^ ) +i-ix-t 



Vp/iq) 



V+\q) 



whose real part is positive. 

The function v : EiU E2 ^^ C is implicitly defined as the only root of the function 



J='{v,q,t) = -z 



1 2V' f I 2V" 

+ v'^ \ + /i — ^ h -y +ivx-t 



Vp/{q) 



V+\q) 
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such that Qm [dtv{t, q)] < 0, with 

El = {{t,q) e m+ X m-^lth, <t<to andO < q < qo{t)} 

and 

E2 = {(t, g) G -K+ X iR+ I to < t < t/ii and qo{t) < q < gi(t)} • 

Up^{x,z,t) is the displacement of the transmitted Ps wave and satisfies: 



^Ps,x(2^'^'0 =-—T / ^^ 

^ Jo 



Up^ Jx,z,t) 



Vl2 r^'^ 



-7r2 , 
^ ./O 



ifthi < t < to and |9m [7(^0, 0)]| < 



^e 
1 

P12 r^^-^^ 



iv{t, q)Tpsiv{t, q)) — {t, q) 



dv 



dq, 



'^Psi'^ii)^<l)'^Psivit, 9))"^(*' 9) 



dg, 



-7r2 , 

vr ./o 



3fie 



P12 r^'^ 



^e 



Up^ (x,z,t) 



^^ Jqo(t) 

P12 Z"*'^*^ 



9^ 
n(.t,q)Tpf{-f{t,q)) — {t,q) 



du 
iv{t,q)Tpsiv{t,q)) — {t,q) 



dq, 
dq 



^ ./o 



3fie 



(i7 



+ ^ r""«e 



^ -'90(4) 

1 



i^Psi^'y'^^ iWpsilit, 9))-^(i, 9) 



di) 

l<'Psi'^it)'<l)'^Psi'^(t, ^^-^ dt *-*' ^^ 






if tQ < t < th2 and l^m [7(^0,0)]! < 



t'rr 



^Ps,x(a^'^'*) 



Pi 2 /'^°^*^ 



-7r2 , 



3fie 



Vs„ 



yx, z,T) „ 



-7r2 , 

vr ./Q 



90(4) 



3f?e 



d'y 
i7(t, q)Tps{-f{t, 9))-^(*> 9) 



(^7 



dg, 



/^psItI*, q))Tps{i{t, <i))-^{t, q) 



dq, 



ifihi < * ^'^'^ 1^"*- [7(^0) 0)]| < — or if tQ < t and |9m [7(^0; 0)]| > 



Vrn 



Kr 



anrf Up^{x,z,t) = e/se. to denotes here the arrival time of the Ps volume wave at 
point (x,0, z), 



thi 



hi 



1 



1 



1 



V- 



1/2 



1 \x\ 

+ 



y^/ ^^ax ' ^max 



(27) 



denotes the arrival time of the Ps head wave at point (x, 0, z) and 

th2 



h^ + z^-hzi^ + ^]+x^ 



h _z_ 

Cl C2 



(28) 
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denotes the time after which there is no longer head wave at point {x,0,z), where 



ci 



V+2 F2 



and C2 



1 1 



1/2 T/2 • 

* Ps max 



The function qq : [to; +cxd] h^ M^ is the reciprocal function of t^ : ]R^ ^^: [to,+oo]; 
where to{q) is the arrival time at point (x,0,z) of the fictitious Ps volume wave, prop- 
agating at a velocity V^{q) in the top layer and at velocity Vpg{q) in the bottom layer. 

The function qi : [ti ; to] ^^ ^^ is defined by 



9l(*) 



\ 



-^ {t + z. 



T/- 2 y2 



T/+2 y2 



V3' 



The function 7 : {{t,q) £ IR^ x IR^ \ t > to(?)} ^^ C is implicitly defined as the only 
root of the function 



^i^,<l,t) = -z 



yps'ii) 



1/2 



+ h 



1 



V+\q) 



1/2 

+ 7^ I + ^73^ ~ t 



whose real part is positive. 

The function v : EiD E2 >-^ C is implicitly defined as the only root of the function 



J^iv,q,t) 



1 



yps'ii) 



\ 1/2 

+ v^\ +h 



1 



V+\q) 



1/2 

+ v'^ \ + ivx — t 



such that 9m [(9it;(t, g)] < 0, with 

El = {{t,q) € IR+ X lR+\th^<t<to andO < q < qo{t)} 

and 

E2 = {{t-,Q) £ IR^ X IR^lto <t <thi and qo{t) < q < qi{t)} . 

Uc{x,z,t) is the displacement of the transmitted S wave and satisfies: 



-7r2 , 



UsJx,Z,t) = -^ / 

^ ^ Jo 



1 fQiit) 



^e 



iv{t, q)Kg{v(t, q))Ts{v{t, Q))-j^{t, q) 



dq, 



9l(*) 



^e 



{v\t,q) + q^)rs{v{t,q))^{t,q) 



dq, 



ifthi <t<tQ and 1 9m [7(^0, 0)] I < 



Kr 
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Us.x(^,Z,t) 



Ug (x,z,t) 



I fQo{t) 






3f?e 



1 r<ii(t) 



n{t, 9)/«5 (7(^1 qWsiiit, Q))-j^iii q) 



di) 
iv{t, q)Kg{v{t, q))Ts{v{t, q))-rr{t, q) 



dt 



dq, 
dq 



-7r2 , 



Qo{t) 



^e 



dj 



1 rii{t) 



■K" 



90 (t) 



(y(t,g) + g^)r5(7(t,g))^(t,c?) 
{v\t,q)+q^)Ts{v{t,q))^{t,q) 



dq 
dq, 



if tQ < t < th^ and l^m [7(^0,0)]! < 



V 

•'max 



1 r 

^ JO 



^e 



_ d'y 

n{i, q)i^s (7(*> (iWs{i{t, 9))"^(*' 9) 



dq, 



Ug {x,z,t) 



I fQoit) 



TT^ 



^e 



{^\t,q) + q')Ts{j{t,q))^{t,q) 



dq, 



if th2 < t and |9m [7(^0, 0)]| < — or if Iq <t and |9m [7(^0; 0)]| > — and 

''max •'max 

Upg{x,z,t) = else, to denotes here the arrival time of the S volume wave at point 
{x,0,z) (we recall in appendix the computation of to), 

h 



thi 



1 



1 



y^ 



1/2 



1 



T/-^ V^ ' K,«x 

V o max * max 



1 \x\ 

+ 



(29) 



denotes the arrival time of the S head-wave at point {x, 0, z) and 



ihn 



h _z_ 

Cl C2 



(30) 



denotes the time after which there is no longer head wave at point (x,0, z), where 



Cl 



V~ 



1/2 



and C2 



Ks ""^max 



The function go : \to; +00] 1-^ IR^ is the reciprocal function of to : SV~ 1— >: [io,+oo], 
where to{q) is the arrival time at point {x,0,z) of the fictitious S volume wave, propa- 
gating at a velocity V^{q) in the top layer and at velocity Vg{q) in the bottom layer (we 
recall in appendix the computation ofto{q)). 

The function qi : \ti ; to] h^ IR^ is defined by 



«i(() 



\ 



-t: \t + z. 



T/-2 1/2 
V Q ' max 



1/+2 1/2 



1/2 • 
'max 
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The function 7 : {{t,q) £ IR^ x ]R^ \t > tQ{q)} 1— > C is implicitly defined as the only 
root of the function 

whose real part is positive. 

The function v : Ei U E2 >-^ C is implicitly defined as the only root of the function 

such that '^m[dtv{t,q)] < 0, with 

El = {{t,q) e m^ X IR+\th^<t<to andO < q < qo{t)} 

and 

E2 = {(*,?) € IR^ X IR^lto <t <th^ and qo{t) < q < qi{t)} . 

3 Proof of the theorem 

To prove the theorem, we use the Cagniard-de Hoop method (see [H EJ [HI [I5l [H] ) , which 
consists of three steps: 

1. We apply a Laplace transform in time, 

»+oo 



u{x,y,z,s)= / u{x,y,z,t)e * dt, 
Jo 

and a Fourier transform in the x and y variables, 

f+00 f+00 
iL{kj:,ky,z,s) = / / u{x,y,s)e'^^''''~^^yy^ dxdy 



00 J —00 



to (J17p in order to obtain an ordinary differential system whose solution Q(kx,ky,z,s) 
can be explicitly computed (§ 13. 1|) : 

2. we apply an inverse Fourier transform in the x and y variables to G (we recall that we 
only need the solution at y = 0: 

g{x,0,z,s) = —^ / g {kx,ky,z,s)e~' '''='' dk^dky. 

And, using tools of complex analysis, we turn the inverse Fourier transform in the x 
variable into the Laplace transform of some function 7i{x^kyjZjt) (§ 13. 2p : 

g{x,0,z,s) = —^ / n{x,ky,z,t)e-'^dtdky; (31) 

47r J_oo Jo 
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3. the last step of the method consists in inverting the order of integration in ()3ip to obtain 



G{x,0,z,s) = -^ 



hq(t) 



7i{x,ky, z,t) dky e * dt. 



'0 \J-q{t) 

Then, using the injectivity of the Laplace transform, we identify Q{x, 0, z, t) to 

1 /■+<?(*) 



47r2 



git) 



/L{3j^ rCy j Z ^ V ] QjfXiij 



(see §E3]). 



3.1 The solution in the Laplace- Fourier plane 

Let us first apply a Laplace transform in time and a Fourier transform in the x and y variables 
to p!7j) to obtain 



2 ~^ k^ + Ky ] p 



V+ 

2 



' +kl + kl\^: 



dz^ V 

a2|>7 



_2 ' 



y>o. 



^_2 ' 'x ' ^y ^ Q. 



0, i€{Pf,Ps,S} y<0, 

y = o, 



(32) 



t ^(p+,cI,-^,$-^,$-) = 
where 6 is the Laplace-Fourier transform of the operator B. 



From the two first equations of ([32]) . we deduce that the solution (p+,(^^ )ie{Pf,Ps,s}) is 
such that 



^ =nnc+^ref' 



1 



P- 



-s\z—h\K~^ 



mc 



^~ = Ti{kj:,ky,s)e v"" 



'k-. k. 



?)). 



" ' ' ^ ^P^q{ = Rikx,ky,s)e 
ie{Pf,Ps,S}, 



' kx ky 



(33) 



where the coefficients R and Tj are computed by using the last equation of ([320 : 

i3(p+f,|.p^,|.p„6^) = -^(p^^^, 0,0,0), 
or, from ( [33]) : 



-shn' 



B [R{k^, ky, s),Tpf{k^, ky, s),Tps{kx, ky, s), rs(A:^, ky, s)] = -B 



SK~' 



+ ( kx ky_ 



-,0,0,0 
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After some calculations that we don't detail here, we obtain that R{kx, ky,s), Tpf{kx, ky,s), 
Tpsikx, ky,s), and Ts{kx, ky,s)) are solution to 



Al^ y 

s ' s 



K[kx, ky, s) 
s^Tpf{kx,ky,s) 

S lps\tix, ky, S) 

s^Ts{kx,ky,s) 



shK+fi^,!^^, 



2SK+ 



kx >^y 



V+' 



l^+ fkx tv 



1 



1 



(34) 



From the definition of the reflection and transmission coefficients we deduce that 

7^ 



li[kx, ky, s) 
s^Tpf{kx,ky,s) 
s 1 ps[kx, ky, s) 
s^Ts{kx,ky,s) 



s ' s 



Tpsi'-f.'-f 



■J- I'ka hi 



-shK+{^) 



(35) 



Finally, we obtain: 



■'^inc ^ref 



P: 



mc 



SV+^K^ 



s ' s 



'^ref 



s V s ' s 



kx ky 

s ' s 



d- = l-r f ^> 1 e 



s ' s 



x '^y \ -SUK. I ^,^)-/iK+(^,^ 



S S 



' kx_ hi 

s ' s 



)). i£{Pf,Ps], 

6- = -Lr^ (^,M e-<^'^^i"''^)''^'^"("'^)). 



(36) 



and 



^ n+ 



U- + U n, 

mc ref 



u- = 1 

mc.a; 



-^P 



«:+(%,% 



p^s 



2^mc' 



u- = siffn(/i — z)- 
mc.z ^ ^ ^ 



p^s 



■p. 



inc' 



ref,x 



kx 



li 



/9"^S 



2^ref' ref, 



^,+ (kx ky 



p^s 



ref 



(37) 



u. 



u. 



-ikxVn^Pf -ikxVu^p^ -iskxKg [ — ) ^^ 



SK 



Pf 



^, ^) Vu^-pf + s.-p^ (^, ^ ) Vu'^-p, + (fe^ + fcP<I> 



s s 



s s 



S)^s 



In the following we only detail the computation of ti^^ p^ = —ikxVu ^p^j since the compu- 
tation of the other terms is very similar. 
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3.2 The Laplace transform of the solution 

We apply an inverse Fourier transform in the x and y variable to uj^ p^ and we set kx = QxS 
and ky = qyS to obtain (we recall that we consider y = 0) 

"OO "/ — oo 

= -^/ / E{qx,qy)dqxdqy, 

J — oo "/ — oo 



-oo J — oo 

with 

S(fe.,gj,) = ig,rp,(g,,gj,)e-<-^'=p^(''-''-)+'^'^+(^-''''))+i^^^). 

Let us now focus on the integral over q^ for a fixed g^^ 

+ 00 /" + 00 _ 

Hfe,g,)dg. = / ig,Tp,(g,,g,)e-<-^'^?=(«-''-)+'^'=^('?-^-»+^^^^)d<7, (38) 

CO J ~oo 

This integral is very similar to the one we have obtained in 2D [8j, therefore, using the same 
method, we have: 

1 



if 17(^2/, *o)| > 



I "Eiqx,qy)dqx = 2 f '^ ?fie (■n{qy,t)Tps{qy,^{qy,t))^^^^^^) e~''dt 

J-oo Jioiqv) V at J 



'^^\i{%M{qy))\ < 



'toiqy) 

1 



'^ma.xyQy ) 

hiqy) 



/ E{qx,qy)dqx = 2 ^e {iv{qy,t)Tps{qy,v{qy,t)) ^^f ^ e'^^dt 

^-oo ^thfe) V ct y 

+ 2 T" 3f?e ('i7(g^,t)Tp,(g„7((z„t))^%^') e-^*di 



where Vmax is greatest fictitious velocity defined by: 



■^iiiax • — "^inaxwi — ''r 



, l + V^iax'g^ 

to is the fictitious arrival time of the Ps volume wave we have defined in the theorem, thiqy) 
the fictitious arrival time of the Ps head wave defined by 



~ ~ 1 1 /I 1 \x\ 

th ■■= th{q) = hJ--T— - TTT^-T - zA-—Y— - „o ,s + 



V+\q) Vl^M ^Vp,\q) Vl^M V^ax(g)' 
Let us recall that the condition \'yiqy,io{qy))\ < — — - is equivalent to 

i^mSuxyQy ) 

l7(0,to)| < TT and \qy\ < gmax, 
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with 



i 



K. 



V ' V^niax 



Moreover, th is bijective from [0 ; gmax] to [to ; t/i('?max)] and we denote its inverse by qh' 
Qhit) 



\ 



1 / 1 

-^ t + z 



1 



Vp^ "ma 



hi 



1 



1 



y+2 y2 



1 

y2" 



Let us also recall that for qy = (?max) the arrival times of the fictitious head and volume waves 
are the same: t/i(c?max) = ^o(9max)- As an illustration, we represent the functions to (the red 
dotted line) and th (the blue solid line) in Fig. [3l 



> 


< t 


io(gmax) = thiqiaax) = th2 


^ 


*0 

th. 


ioiqy) ....■■■■^^^ 




---^ th{qy) 


1 >► 



Figure 3: Functions to (red, dotted line) and th (blue, solid line) 

We then deduce that 
. if |7(0,to)| >^ 



4h(,.,,,) 

if|7(0,to)|<^ 



dq^dqy = 2 \ \ ^e [ij{qy,t)Tps{qy,j{qy,t)) 

MJioiQy) V dt 



% 



^^(^^'*)^e--dtdg,; 



I E{q^,qy)dqxdqy = 2 / ^e {iv{qy,t)Tps{qy,v{qy,t)) ^^ ] e'"" dtdqy 

^IR ^-gmax Jihiqy) V Ct / 

+ 2/ f"^ ^e(n{qy,t)rps{qy,j{qy,t))^^:^^^)e-''dtdqy. 
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3.3 Inversion of the integrals 

The key point of the method is the inversion of the integral with respect to qy with the integral 
with respect to t. For the volume wave we have (see Figs. [Hand [5]), after having remark that 
the integrand is even with respect to qy: 



/ + 00 /" + 00 / 

/ ?lieii-f{qy,t)Tps{qy,-fiqy,t)) 
-oo Jin(n-n) \ 



to(qy) 
+00 rqo(i) 



^^^"^y^'^^e-^dtdqy 



Ke i-f{qy,t)Tps{qy,-f{qy,t)) 



dt 



^^]e--dqydt; 



'to JO 
and for the head wave (see Figs. El and [7]): 

Ke iv{qy,t)Tps{qy,v{qy,t)) ^^J' ^ e-'' dtdqy 



9max Jtfi{qy) 

to rqh{t) 



th. JO 



3f?e iv{qy,t)Tps{qy,v{qy,t)) 



dt 



dv{qy,t) 
dt 



e dqydt 



+ 2['' r ^e(iv{qy,t)Tps{qy,viqy,t))^^^%^)e-'' dqydt. 

J to Jqo(t) \ 



dt 





Figure 4: Integration first over 
t for the volume wave 



qy qy 

then over Figure 5: Integration first over t then over qy 
for the volume wave 




to 
t/ii 



io(9max) — thiQra&yi) — ih2 




Qy Qy 

Figure 6: Integration first over qy then over Figure 7: Integration first over t then over qy 

t for the head wave for the head wave 
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We thus have: 

Loo 

-St 



^sa;,Ps(^'0'^''5) = / u^^p^{x,0,z,t)e ' dt 
Jo 

and we conclude by using the injectivity of the Laplace transform. 

4 Numerical illustration 

To illustrate our results, we have computed the green function and the analytical solution to 
the following problem: we consider an acoustic layer with a density /o"*" = 1020 kg/m^ and a 
celerity V~^ = 1500 m/s on top of a poroelastic layer whose characteristic coefficients are: 

• the solid density p~ = 2500 kg/m^; 

• the fluid density p7 = 1020 kg/m^; 

• the porosity cj}~ = 0.4; 

• the tortuosity a~ = 2; 

• the solid bulk modulus K' = 16.0554 GPa; 

• the fluid bulk modulus KJ = 2.295 GPa; 

• the frame bulk modulus K/^ = 10 GPa; 

• the frame shear modulus ;U^ = 9.63342 GPa; 

so that the celerity of the waves in the poroelastic medium are: 

• for the fast P wave, Vpr = 3677 m/s 

• for the slow P wave, Vf^ = 1060 m/s 

• for the V' wave, Vg = 2378 m/s. 

The source is located in the acoustic layer, at 500 m from the interface. It is a point source 
in space and a fifth derivative of a Gaussian of dominant frequency /o = 15 Hz: 



fit) = 2 ,2 
JO 



vr- „ ..vr^ /. 1\^ .vr^ /. 1\' "t^I*-^ 



3 + 12-2 t-- +4-jt 

/o V foj /o V /o 



7JV /o 



We compute the solution at two receivers, the first one is in the acoustic layer, at 533 m from 
the interface; the first one is in the poroelastic layer, at 533 m from the interface; both are 
located on a vertical line at 400 m from the source (see Fig. [8|). To compute the integrals 
over q and the convolution with the source function, we used a classical mid-point quadrature 
formula. 

We represent the z component of the green function associated to the displacement from 
i = 0tot = 1.2son Fig [9] and the displacement in Fig. [lOl The left picture represents 
the solution at receiver 1 while the right picture represents the solution at receiver 2. As all 
the types of waves are computed independently, it is easy to distinguish all of them, as it is 
indicated in the figures, solution. 
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Figure 8: Configuration of the experiment 
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Figure 9: The z component of the green function associated to the displacement at receiver 
1 (left picture) and 2 (right picture). 
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Figure 10: The z component of the displacement at receiver 1 (left picture) and 2 (right 
picture). 
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5 Conclusion 

In this paper we have provided the complete solution (reflected and transmitted wave) of the 
propagation of wave in a stratified 3D medium composed of an acoustic and a poroelastic 
layer. In a forthcoming paper we will extend the method to the propagation of waves in 
bilayered poroelastic medium in three dimensions. 

A Definition of the fictitious and real arrival times of the vol- 
ume waves. 

We detail in this section the computation of the fictitious and real arrival times of the trans- 
mitted Ps wave at point {x,0,z). For a given qy G IR, we first determine fastest path of the 
wave from the source to the point {x,0,z), travelling at a velocity V~^{qy) in the upper layer 
and at a velocity Vp^i^Qy) in the bottom layer: we search a point ,^0 on the interface between 
the two media which minimizes the function 



^^' v+(q,) + VpM 

(see Fig. [TTI) . This leads us to find ^0 such that 

i'(^o) = ^° + ^°~^ = 0. (39) 

From a numerical point of view, the solution of this equation is done by computing the roots 
of the following fourth degree polynomial 



1 1 \ „. „ / 1 1 



V+'iqy) V-p^\qy)l \Vps\%) ^^'i%) , 



X^ + 2x\ . . X^ 



x^ + 2;^ x"^ + h? \ ,,9 xh? ^^ x^h? 

X^H ^ X + 



\V+\qy) Vp^\qy)J Vpfiqy) Vp^^qy)' 

^0 is thus the only real root of this polynomial located between and x which is also solution 
of ([39|) . Once ^0 is computed, we can define 



toiqy) = -^J^^ + V^.(,,) ^"' '' = ''^'^ 

Let us remark that 

Property A.l. Since the fictitious velocities are smaller than the real one, the fictitious 
arrival times are greater than the real one. Moreover, since the fictitious velocities are even 
functions decreasing on ]R^ , Iq is an even function, increasing on ]R^ . 

Corollaire A.l. The function to is bijective from IR^ to IR^ . 
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Figure 11: Path of the transmitted i wave 
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